Active Contours using Parameteric Curves
This tour explores image segmentation using parametric active contours.
Contents
Installing toolboxes and setting up the path.
You need to download the following files: signal toolbox, general toolbox and graph toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_signal, toolbox_general and toolbox_graph in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(p,path); % scilab users must *not* execute this
Then you can add the toolboxes to the path.
getd('toolbox_signal/'); getd('toolbox_general/'); getd('toolbox_graph/');
Parameteric Curves
In this tours, the active contours are represented using parametric curve. It is implemented using a picewise linear curve, which is just a vector of complex numbers.
We compute an initial curve by performing subdivision of a control polygon.
Initial polygon.
c0 = [0.78 0.14 0.42 0.18 0.32 0.16 0.75 0.83 0.57 0.68 0.46 0.40 0.72 0.79 0.91 0.90]' + ...
1i* [0.87 0.82 0.75 0.63 0.34 0.17 0.08 0.46 0.50 0.25 0.27 0.57 0.73 0.57 0.75 0.79]';
Interpolate the curve.
p = 256;
c1 = c0;
c1(end+1) = c1(1);
d = abs(c1(1:end-1)-c1(2:end)); d = [0;cumsum(d)];
c1 = interp1(d/d(end),c1,(0:p-1)'/p, 'linear');
Display the initial curve.
clf; h = plot(c1([1:end 1]), 'k'); set(h, 'LineWidth', 2); axis('tight'); axis('off');
Compute the normal to the curve. This is obtained by rotating by pi/2 the tangent.
tgt = c1([2:end 1]) - c1([end 1:end-1]); normal = -(1i*tgt)./abs(tgt);
Move the curve in the normal direction.
h = .03; c2 = c1 + normal*h; c3 = c1 - normal*h;
Display the curves.
clf; hold on; h = plot(c1([1:end 1]), 'k'); set(h, 'LineWidth', 2); h = plot(c2([1:end 1]), 'r--'); set(h, 'LineWidth', 2); h = plot(c3([1:end 1]), 'b--'); set(h, 'LineWidth', 2); axis('tight'); axis('off');
Evolution by Mean Curvature
The simplest evolution is the mean curvature evolution, that is the evolution that reduces the length of the curve.
Time step for the evolution.
dt = 0.001;
Initialize the curve.
c = c1;
Compute the tangent using forward derivatives.
e1 = c([2:end 1]) - c; e1 = e1 ./ abs(e1);
Compute the normal time curvature using backward derivatives.
e2 = e1 - e1([end 1:end-1]);
Evolution of the curve.
c = c + dt * e2;
To stabilize the evolution, it is important to re-sample the curve so that it is unit-speed parametrized. You do not need to do it every time step though (to speed up).
c(end+1) = c(1); d = abs(c(1:end-1)-c(2:end)); d = [0;cumsum(d)]; c = interp1(d/d(end),c,(0:p-1)'/p);
Exercice 1: (the solution is exo1.m) Perform the curve evolution, for a time of Tmax=2. You need to resample it a few times.
exo1;
Geodesic Active Contours
Geodesic active contours minimize a weighted length, where the weight is measured according to some image.
Create a synthetic image with some dots.
n = 200; nbumps = 40; theta = rand(nbumps,1)*2*pi; r = .6*n/2; a = [.62*n .6*n]; x = round( a(1) + r*cos(theta) ); y = round( a(2) + r*sin(theta) ); W = zeros(n); W( x + (y-1)*n ) = 1; W = perform_blurring(W,10); W = rescale( -min(W,.05), .3,1);
Display the metric.
clf; imageplot(W);
Pre-compute the derivatives.
options.order = 2; G = grad(W, options); G = G(:,:,1) + 1i*G(:,:,2);
Initialize the curve with a circle.
r = .98*n/2;
p = 128; % number of points on the curve
theta = linspace(0,2*pi,p+1)'; theta(end) = [];
c0 = n/2*(1+1i) + r*(cos(theta) + 1i*sin(theta));
c = c0;
For this experiment, the time step should be larger (because the curve is in [1,n] x [1,n])
dt = .8;
Display the curve on the back ground;
lw = 2; clf; hold on; imageplot(W); h = plot(imag(c([1:end 1])),real(c([1:end 1])), 'r'); set(h, 'LineWidth', lw); axis('ij');
Evaluate the gradient.
g = interp2(1:n,1:n, G, imag(c), real(c));
Evaluate the potential.
w = interp2(1:n,1:n, W, imag(c), real(c));
Compute the tangent using forward derivatives.
e1 = c([2:end 1]) - c; d1 = abs(e1); e1 = e1 ./ d1;
Compute the normal time curvature using backward derivatives.
e2 = w.*e1 - w([end 1:end-1]).*e1([end 1:end-1]);
Evolution of the curve.
c = c + dt * ( - g .* d1 + e2 );
Exercice 2: (the solution is exo2.m) Perform the curve evolution, for a time of Tmax=1150. Remeber to re-sample the curve several time during the evolution.
exo2;
Medical Image Segmentation
One can use a gradient-based metric to perform edge detection in medical images.
Load an image.
n = 256;
M = rescale( sum(load_image('cortex', n), 3 ) );
Display.
clf; imageplot(M);
Exercice 3: (the solution is exo3.m) Compute an edge attracting criterion W, that is small in area of strong gradient. You can use, among other, the function grad, perform_blurring, and threshold too large gradients.
exo3;
Exercice 4: (the solution is exo4.m) Create an initial circle c0 of p=128 points.
exo4;
Exercice 5: (the solution is exo5.m) Perform the curve evolution, for a time of Tmax=1150. Remeber to re-sample the curve several time during the evolution. Try with different dynamics for the W scaling so that you capture the contour you are interested in.
exo5;
Evolution of a Non-closed Curve
It is possible to perform the evolution of a non-closed curve by adding boundary constraint (starting and ending points).
In this case, the algorithm find a local minimizer of the geodesic distance between the two points.
Note that a much more efficient way to solve this problem is to use the Fast Marching algorithm to find the global minimizer of the geodesic length.
Load an image.
n = 256;
M = rescale( sum(load_image('cortex', n), 3 ) );
M = M(46:105,61:120);
n = size(M,1);
Display.
imageplot(M);
Exercice 6: (the solution is exo6.m) Compute an edge attracting criterion W, that is small in area of strong gradient. You can use, among other, the function grad, perform_blurring, and threshold too large gradients.
exo6;
Start and end points.
pstart = 4 + 55i; pend = 53 + 4i;
Initial curve.
p = 128; t = linspace(0,1,p)'; c0 = t*pend + (1-t)*pstart;
Initialize the evolution.
c = c0;
Display.
clf; hold on; imageplot(W); h = plot(imag(c([1:end])),real(c([1:end])), 'r'); set(h, 'LineWidth', 2); h = plot(imag(c([1 end])),real(c([1 end])), 'b.'); set(h, 'MarkerSize', 30); axis('ij');
Exercice 7: (the solution is exo7.m) Perform the curve evolution, for a time of Tmax=1000. Be careful to impose the boundary conditions at each step. Remeber to re-sample the curve several time during the evolution.
exo7;
